################################################################################
### Elite Responses to Ethnic Diversity and Interethnic Contact              ###
### Analysis File                                                            ###
### William O'Brochta                                                        ###
### Louisiana Tech University                                                ###
################################################################################

library(plyr)
library(dplyr)
library(lme4)
library(brms)
library(lmerTest)
library(lmtest)
library(MASS)
library(sandwich)
library(ordinal)
library(multiwayvcov)
library(texreg)
library(mediation)
library(psych)


#Survey3.2 with affect
#load("survey3.2.RData")
source("InterethnicContactPrep.R")

#Descriptive statistics
rm(gender)
attach(survey3.2)
data_summary<-as.data.frame(cbind(cmte_hh_index, cmte_pct_notforward,
                                  contact_avg, contactf_avg,
                                  cmte_spend_wisely, cmte_opinions,
                                  cmte_validconcerns, cmte_onegroup,
                                  cmte_castetrust, talk_comfort, neighbor_comfort,
                                  cmte_interaction_pos01, enthusiastic,
                                  angry, hopeful, resentful, donation_caste,
                                  age, gender, year_elected, caste_reservation, Brahmin,
                                  OF, OBC, SC, ST, Other_Religion, BJP, INC,
                                  edu_high, years_6over, socialmedia_active, 
                                  times_called_3over))
  
detach(survey3.2)
data_summary2<-as.data.frame(matrix(nrow=33, ncol=6))
for(i in 1:33){
  data_summary2[i,1]<-names(data_summary)[i]
  data_summary2[i,2]<-min(data_summary[,i], na.rm=T)
  data_summary2[i,3]<-max(data_summary[,i], na.rm=T)
  data_summary2[i,4]<-sqrt(var(data_summary[,i], na.rm=T))
  data_summary2[i,5]<-mean(data_summary[,i], na.rm=T)
  data_summary2[i,6]<-median(data_summary[,i], na.rm=T)
}
colnames(data_summary2)<-c("Variable", "Min", "Max", "SD", "Mean", "Median")

###Table SI.1.1
print(xtable(data_summary2, type='latex', digits=2), file='data_summary2.tex', include.rownames=FALSE)

#Descriptives at committee level
survey3.2$cmte_surveynum<-1
committee_summary<-as.data.frame(cbind(survey3.2$committee_asked_about,
                                       survey3.2$cmte_size,
                                       data_summary))
names(committee_summary)[1]<-"committee_asked_about"
names(committee_summary)[2]<-"cmte_size"

committee_summary2<-committee_summary%>%
  group_by(committee_asked_about)%>%
  summarise_all(mean, na.rm=T)

committee_summary.1<-as.data.frame(cbind(survey3.2$committee_asked_about,
                                         survey3.2$cmte_surveynum))
names(committee_summary.1)[1]<-"committee_asked_about"
names(committee_summary.1)[2]<-"cmte_surveynum"
committee_summary.1$cmte_surveynum<-as.numeric(as.character(committee_summary.1$cmte_surveynum))

committee_summary2.1<-committee_summary.1%>%
  group_by(committee_asked_about)%>%
  summarise_all(sum, na.rm=T)

committee_summary2$cmte_surveynum<-committee_summary2.1$cmte_surveynum
committee_summary2<-as.data.frame(committee_summary2)

committee_summary3<-as.data.frame(matrix(nrow=35, ncol=6))
for(i in 1:35){
  committee_summary3[i,1]<-names(committee_summary2)[i+1]
  committee_summary3[i,2]<-min(committee_summary2[,i+1], na.rm=T)
  committee_summary3[i,3]<-max(committee_summary2[,i+1], na.rm=T)
  committee_summary3[i,4]<-sqrt(var(committee_summary2[,i+1], na.rm=T))
  committee_summary3[i,5]<-mean(committee_summary2[,i+1], na.rm=T)
  committee_summary3[i,6]<-median(committee_summary2[,i+1], na.rm=T)
}
colnames(committee_summary3)<-c("Variable", "Min", "Max", "SD", "Mean", "Median")

#print(xtable(committee_summary3, type='latex', digits=2), file='committee_summary.tex', include.rownames=FALSE)

#Range of diversity on committees in corporation with more than one committee
committee_summary3.1<-committee_summary2
committee_summary3.1$cmte_totalnum<-1
committee_summary3.1$corporation<-gsub("_.*", "", committee_summary3.1$committee_asked_about)
committee_summary3.1<-dplyr::select(committee_summary3.1, c("corporation", "cmte_totalnum", "cmte_hh_index"))
committee_summary3.2<-as.data.frame(committee_summary3.1%>%
  group_by(corporation)%>%
  summarise_each(list(min=min,max=max,sum=sum)))
committee_summary3.2<-committee_summary3.2[committee_summary3.2$cmte_totalnum_sum>1,]
committee_summary3.2$range<-committee_summary3.2$cmte_hh_index_max-committee_summary3.2$cmte_hh_index_min
mean(committee_summary3.2$range)


#Caste Coding Comparison

#Plot for cmte_hh_index
mod<-lm(cmte_hh_index_sr~cmte_hh_index, survey3.2)
newx <- seq(min(survey3.2$cmte_hh_index), max(survey3.2$cmte_hh_index), length.out=100)
preds <- predict(mod, newdata = data.frame(cmte_hh_index=newx), 
                 interval = 'confidence')

plot(jitter(cmte_hh_index_sr, amount=.015)~jitter(cmte_hh_index, amount=.015), survey3.2, lwd=2, 
     xlab="Coded Diversity", ylab="Self-Reported Diversity",
     xlim=c(0,0.8), ylim=c(0,0.8), main="Coded and Self-Reported Diversity",
     cex.lab=1.3, cex.axis=1.3, cex.main=1.3)
abline(mod, lwd=3)
lines(newx, preds[ ,3], lty = 'dashed', col = 'red', lwd=3)
lines(newx, preds[ ,2], lty = 'dashed', col = 'red', lwd=3)

cor(survey3.2$cmte_hh_index, survey3.2$cmte_hh_index_sr, use="complete.obs")


#Plot for contact_avg
mod1<-lm(contact_avg_sr~contact_avg, survey3.2)
newx1 <- seq(0,1, length.out=100)
preds1 <- predict(mod1, newdata = data.frame(contact_avg=newx1), 
                 interval = 'confidence')

plot(jitter(contact_avg_sr, amount=.015)~jitter(contact_avg, amount=.015), survey3.2, lwd=2, 
     xlab="Coded Contact", ylab="Self-Reported Contact",
     xlim=c(0,1), ylim=c(0,1), main="Coded and Self-Reported Contact",
     cex.lab=1.3, cex.axis=1.3, cex.main=1.3)
abline(mod1, lwd=3)
lines(newx1, preds1[ ,3], lty = 'dashed', col = 'red', lwd=3)
lines(newx1, preds1[ ,2], lty = 'dashed', col = 'red', lwd=3)

cor(survey3.2$contact_avg, survey3.2$contact_avg_sr, use="complete.obs")


###Check random assignment to treatment
library(nnet)
library(stargazer)

#Wald Test
model0<-lm(treated~female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                     caste_reservation+corporation, survey3.2)
waldtest(model0)

#Table SI.3.1
model0.1<-multinom(treated~
                     female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                     caste_reservation+corporation, survey3.2)
z.1 <- summary(model0.1)$coefficients/summary(model0.1)$standard.errors
(p.1 <- (1 - pnorm(abs(z.1), 0, 1)) * 2)
stargazer(model0.1)


#Table SI.3.2
(model0.2<-lm(treated~female, data=survey3.2))
mean(survey3.2[survey3.2$treated==1,]$female, na.rm=T)
mean(survey3.2[survey3.2$treated==0,]$female, na.rm=T)
(model0.3<-lm(treated~age, data=survey3.2))
mean(survey3.2[survey3.2$treated==1,]$age, na.rm=T)
mean(survey3.2[survey3.2$treated==0,]$age, na.rm=T)
(model0.4<-lm(treated~socialmedia_active, data=survey3.2))
mean(survey3.2[survey3.2$treated==1,]$socialmedia_active, na.rm=T)
mean(survey3.2[survey3.2$treated==0,]$socialmedia_active, na.rm=T)
(model0.5<-lm(treated~edu_high, data=survey3.2))
mean(survey3.2[survey3.2$treated==1,]$edu_high, na.rm=T)
mean(survey3.2[survey3.2$treated==0,]$edu_high, na.rm=T)
(model0.6<-lm(treated~years_6over, data=survey3.2))
mean(survey3.2[survey3.2$treated==1,]$years_6over, na.rm=T)
mean(survey3.2[survey3.2$treated==0,]$years_6over, na.rm=T)
(model0.7<-lm(treated~BJP, data=survey3.2))
mean(survey3.2[survey3.2$treated==1,]$BJP, na.rm=T)
mean(survey3.2[survey3.2$treated==0,]$BJP, na.rm=T)
(model0.8<-lm(treated~INC, data=survey3.2))
mean(survey3.2[survey3.2$treated==1,]$INC, na.rm=T)
mean(survey3.2[survey3.2$treated==0,]$INC, na.rm=T)
(model0.9<-lm(treated~times_called_3over, data=survey3.2))
mean(survey3.2[survey3.2$treated==1,]$times_called_3over, na.rm=T)
mean(survey3.2[survey3.2$treated==0,]$times_called_3over, na.rm=T)
(model0.10<-lm(treated~caste_reservation, data=survey3.2))
mean(survey3.2[survey3.2$treated==1,]$caste_reservation, na.rm=T)
mean(survey3.2[survey3.2$treated==0,]$caste_reservation, na.rm=T)

(model0.11<-lm(treated~Brahmin, data=survey3.2))
mean(survey3.2[survey3.2$treated==1,]$Brahmin, na.rm=T)
mean(survey3.2[survey3.2$treated==0,]$Brahmin, na.rm=T)
(model0.12<-lm(treated~OF, data=survey3.2))
mean(survey3.2[survey3.2$treated==1,]$OF, na.rm=T)
mean(survey3.2[survey3.2$treated==0,]$OF, na.rm=T)
(model0.13<-lm(treated~SC, data=survey3.2))
mean(survey3.2[survey3.2$treated==1,]$SC, na.rm=T)
mean(survey3.2[survey3.2$treated==0,]$SC, na.rm=T)
(model0.14<-lm(treated~ST, data=survey3.2))
mean(survey3.2[survey3.2$treated==1,]$ST, na.rm=T)
mean(survey3.2[survey3.2$treated==0,]$ST, na.rm=T)
(model0.15<-lm(treated~OBC, data=survey3.2))
mean(survey3.2[survey3.2$treated==1,]$OBC, na.rm=T)
mean(survey3.2[survey3.2$treated==0,]$OBC, na.rm=T)
(model0.16<-lm(treated~Other_Religion, data=survey3.2))
mean(survey3.2[survey3.2$treated==1,]$Other_Religion, na.rm=T)
mean(survey3.2[survey3.2$treated==0,]$Other_Religion, na.rm=T)






###Analysis

#Table SI.6.1
######## Add in treatment separately
model11<-lm(affect_mr_b2_p~treated+cmte_hh_index+contact_avg+
              female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
              caste_reservation+corporation, survey3.2)
(coef11<-coeftest(model11, vcov=function(x) cluster.vcov(x, survey3.2$corporation)))
marginal<-as.data.frame(summary(margins::margins(model11, 
                                                 variables=c("treated", "cmte_hh_index", "contact_avg"), 
                                                 vcov=vcovCL(model11, cluster=survey3.2$corporation))))
marginal$dv<-"Positive"


model11.1<-glmer(affect_mr_b2_pf~treated+cmte_hh_index+contact_avg+
                   female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                   caste_reservation+(1|corporation), family="binomial",
                 data=survey3.2, control=glmerControl(optimizer = "bobyqa"))
summary(model11.1)




model11.2<-lm(affect_mr_b2_u~treated+cmte_hh_index+contact_avg+
                female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                caste_reservation+corporation, survey3.2)
(coef11.2<-coeftest(model11.2, vcov=function(x) cluster.vcov(x, survey3.2$corporation)))
marginal_temp<-as.data.frame(summary(margins::margins(model11.2, 
                                                      variables=c("treated", "cmte_hh_index", "contact_avg"), 
                                                      vcov=vcovCL(model11.2, cluster=survey3.2$corporation))))
marginal_temp$dv<-"Unpleasant"
marginal<-rbind(marginal, marginal_temp)



model11.3<-glmer(affect_mr_b2_uf~treated+cmte_hh_index+contact_avg+
                   female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                   caste_reservation+(1|corporation),family="binomial",
                 data=survey3.2, control=glmerControl(optimizer = "bobyqa"))
summary(model11.3)


model11.4<-lm(affect_mr_b2_m~treated+cmte_hh_index+contact_avg+
                female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                caste_reservation+corporation, survey3.2)
(coef11.4<-coeftest(model11.4, vcov=function(x) cluster.vcov(x, survey3.2$corporation)))
marginal_temp<-as.data.frame(summary(margins::margins(model11.4, 
                                                      variables=c("treated", "cmte_hh_index", "contact_avg"), 
                                                      vcov=vcovCL(model11.4, cluster=survey3.2$corporation))))
marginal_temp$dv<-"Mixed"
marginal<-rbind(marginal, marginal_temp)



model11.5<-glmer(affect_mr_b2_mf~treated+cmte_hh_index+contact_avg+
                   female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                   caste_reservation+(1|corporation),family="binomial",
                 data=survey3.2, control=glmerControl(optimizer = "bobyqa"))
summary(model11.5)

model11.6<-lm(affect_mr_b2_w~treated+cmte_hh_index+contact_avg+
                female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                caste_reservation+corporation, survey3.2)
(coef11.6<-coeftest(model11.6, vcov=function(x) cluster.vcov(x, survey3.2$corporation)))
marginal_temp<-as.data.frame(summary(margins::margins(model11.6, 
                                                      variables=c("treated", "cmte_hh_index", "contact_avg"), 
                                                      vcov=vcovCL(model11.6, cluster=survey3.2$corporation))))
marginal_temp$dv<-"Weak"
marginal<-rbind(marginal, marginal_temp)



model11.7<-glmer(affect_mr_b2_wf~treated+cmte_hh_index+contact_avg+
                   female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                   caste_reservation+(1|corporation),family="binomial",
                 data=survey3.2, control=glmerControl(optimizer = "bobyqa"))
summary(model11.7)




model11.8<-lm(cmte_spend_wisely_n~treated+cmte_hh_index+contact_avg+
                female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                caste_reservation+corporation, survey3.2)
(coef11.8<-coeftest(model11.8, vcov=function(x) cluster.vcov(x, survey3.2$corporation)))
marginal_temp<-as.data.frame(summary(margins::margins(model11.8, 
                                                      variables=c("treated", "cmte_hh_index", "contact_avg"), 
                                                      vcov=vcovCL(model11.8, cluster=survey3.2$corporation))))
marginal_temp$dv<-"SpendWisely"
marginal<-rbind(marginal, marginal_temp)



model11.9<-ordinal::clmm(cmte_spend_wiselyf~treated+cmte_hh_index+contact_avg+
                           female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                           caste_reservation+(1|corporation),
                         data=survey3.2)
summary(model11.9)


model11.10<-lm(cmte_opinions_n~treated+cmte_hh_index+contact_avg+
                 female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                 caste_reservation+corporation, survey3.2)
(coef11.10<-coeftest(model11.10, vcov=function(x) cluster.vcov(x, survey3.2$corporation)))
marginal_temp<-as.data.frame(summary(margins::margins(model11.10, 
                                                      variables=c("treated", "cmte_hh_index", "contact_avg"), 
                                                      vcov=vcovCL(model11.10, cluster=survey3.2$corporation))))
marginal_temp$dv<-"Opinions"
marginal<-rbind(marginal, marginal_temp)


model11.11<-ordinal::clmm(cmte_opinionsf~treated+cmte_hh_index+contact_avg+
                            female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                            caste_reservation+(1|corporation),
                          data=survey3.2)
summary(model11.11)




model11.12<-lm(cmte_validconcerns_n~treated+cmte_hh_index+contact_avg+
                 female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                 caste_reservation+corporation, survey3.2)
(coef11.12<-coeftest(model11.12, vcov=function(x) cluster.vcov(x, survey3.2$corporation)))
marginal_temp<-as.data.frame(summary(margins::margins(model11.12, 
                                                      variables=c("treated", "cmte_hh_index", "contact_avg"), 
                                                      vcov=vcovCL(model11.12, cluster=survey3.2$corporation))))
marginal_temp$dv<-"ValidConcerns"
marginal<-rbind(marginal, marginal_temp)


model11.13<-ordinal::clmm(cmte_validconcernsf~treated+cmte_hh_index+contact_avg+
                            female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+caste_reservation+
                            (1|corporation),
                          data=survey3.2)
summary(model11.13)



model11.14<-lm(cmte_onegroup_n~treated+cmte_hh_index+contact_avg+
                 female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                 caste_reservation+corporation, survey3.2)
(coef11.14<-coeftest(model11, vcov=function(x) cluster.vcov(x, survey3.2$corporation)))
marginal_temp<-as.data.frame(summary(margins::margins(model11.14, 
                                                      variables=c("treated", "cmte_hh_index", "contact_avg"), 
                                                      vcov=vcovCL(model11.14, cluster=survey3.2$corporation))))
marginal_temp$dv<-"OneGroup"
marginal<-rbind(marginal, marginal_temp)


model11.15<-ordinal::clmm(cmte_onegroupf~treated+cmte_hh_index+contact_avg+
                            female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                            caste_reservation+(1|corporation),
                          data=survey3.2)
summary(model11.15)




model11.16<-lm(cmte_castetrust_n~treated+cmte_hh_index+contact_avg+
                 female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                 caste_reservation+corporation, survey3.2)
(coef11.16<-coeftest(model11.16, vcov=function(x) cluster.vcov(x, survey3.2$corporation)))
marginal_temp<-as.data.frame(summary(margins::margins(model11.16, 
                                                      variables=c("treated", "cmte_hh_index", "contact_avg"), 
                                                      vcov=vcovCL(model11.16, cluster=survey3.2$corporation))))
marginal_temp$dv<-"CasteTrust"
marginal<-rbind(marginal, marginal_temp)


model11.17<-ordinal::clmm(cmte_castetrustf~treated+cmte_hh_index+contact_avg+
                            female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                            caste_reservation+(1|corporation),
                          data=survey3.2)
summary(model11.17)




model11.18<-lm(neighbor_comfort_n~treated+cmte_hh_index+contact_avg+
                 female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                 caste_reservation+corporation, survey3.2)
(coef11.18<-coeftest(model11.18, vcov=function(x) cluster.vcov(x, survey3.2$corporation)))
marginal_temp<-as.data.frame(summary(margins::margins(model11.18, 
                                                      variables=c("treated", "cmte_hh_index", "contact_avg"), 
                                                      vcov=vcovCL(model11.18, cluster=survey3.2$corporation))))
marginal_temp$dv<-"Neighbor"
marginal<-rbind(marginal, marginal_temp)


model11.19<-ordinal::clmm(neighbor_comfortf~treated+cmte_hh_index+contact_avg+
                            female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                            caste_reservation+(1|corporation),
                          data=survey3.2)
summary(model11.19)


model11.20<-lm(talk_comfort_n~treated+cmte_hh_index+contact_avg+
                 female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                 caste_reservation+corporation, survey3.2)
(coef11.20<-coeftest(model11.20, vcov=function(x) cluster.vcov(x, survey3.2$corporation)))
marginal_temp<-as.data.frame(summary(margins::margins(model11.20, 
                                                      variables=c("treated", "cmte_hh_index", "contact_avg"), 
                                                      vcov=vcovCL(model11.20, cluster=survey3.2$corporation))))
marginal_temp$dv<-"Talk"
marginal<-rbind(marginal, marginal_temp)


model11.21<-ordinal::clmm(talk_comfortf~treated+cmte_hh_index+contact_avg+
                            female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                            caste_reservation+(1|corporation),
                          data=survey3.2)
summary(model11.21)


model11.22<-lm(donation_caste~treated+cmte_hh_index+contact_avg+
                 female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                 caste_reservation+corporation, survey3.2)
(coef11.22<-coeftest(model11.22, vcov=function(x) cluster.vcov(x, survey3.2$corporation)))
marginal_temp<-as.data.frame(summary(margins::margins(model11.22, 
                                                      variables=c("treated", "cmte_hh_index", "contact_avg"), 
                                                      vcov=vcovCL(model11.22, cluster=survey3.2$corporation))))
marginal_temp$dv<-"Donation"
marginal<-rbind(marginal, marginal_temp)


model11.23<-glmer(donation_castef~treated+cmte_hh_index+contact_avg+
                    female+age+BJP+INC+
                    (1|corporation), family="binomial",
                  data=survey3.2, control=glmerControl(optimizer = "bobyqa"))
summary(model11.23)



#Linear Models SI.6.1 (main results)
library(stargazer)
stargazer(model11, model11.2, model11.4, model11.6, model11.8, model11.10, model11.12,
          model11.14, model11.16, model11.18, model11.20, model11.22,
          se=list(coef11[,2], coef11.2[,2], coef11.4[,2], coef11.6[,2], coef11.8[,2], coef11.10[,2],
                  coef11.12[,2], coef11.14[,2], coef11.16[,2], coef11.18[,2], coef11.20[,2], coef11.22[,2]))


#Glmer and CLMM models SI.6.2
texreg(list(model11.1, model11.3, model11.5, model11.7, model11.9, model11.11, model11.13,
            model11.15, model11.17, model11.19, model11.21, model11.23), stars = c(0.01, 0.05, 0.1))

marginal7<-marginal




#Table SI.6.5
#Now with state FE
model12<-lm(affect_mr_b2_p~treated+cmte_hh_index+contact_avg+
              female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
              caste_reservation+state, survey3.2)
(coef12<-coeftest(model12, vcov=function(x) cluster.vcov(x, survey3.2$state)))


model12.2<-lm(affect_mr_b2_u~treated+cmte_hh_index+contact_avg+
                female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                caste_reservation+state, survey3.2)
(coef12.2<-coeftest(model12.2, vcov=function(x) cluster.vcov(x, survey3.2$state)))

model12.4<-lm(affect_mr_b2_m~treated+cmte_hh_index+contact_avg+
                female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                caste_reservation+state, survey3.2)
(coef12.4<-coeftest(model12.4, vcov=function(x) cluster.vcov(x, survey3.2$state)))

model12.6<-lm(affect_mr_b2_w~treated+cmte_hh_index+contact_avg+
                female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                caste_reservation+state, survey3.2)
(coef12.6<-coeftest(model12.6, vcov=function(x) cluster.vcov(x, survey3.2$state)))

model12.8<-lm(cmte_spend_wisely_n~treated+cmte_hh_index+contact_avg+
                female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                caste_reservation+state, survey3.2)
(coef12.8<-coeftest(model12.8, vcov=function(x) cluster.vcov(x, survey3.2$state)))


model12.10<-lm(cmte_opinions_n~treated+cmte_hh_index+contact_avg+
                 female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                 caste_reservation+state, survey3.2)
(coef12.10<-coeftest(model12.10, vcov=function(x) cluster.vcov(x, survey3.2$state)))


model12.12<-lm(cmte_validconcerns_n~treated+cmte_hh_index+contact_avg+
                 female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                 caste_reservation+state, survey3.2)
(coef12.12<-coeftest(model12.12, vcov=function(x) cluster.vcov(x, survey3.2$state)))

model12.14<-lm(cmte_onegroup_n~treated+cmte_hh_index+contact_avg+
                 female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                 caste_reservation+state, survey3.2)
(coef12.14<-coeftest(model12, vcov=function(x) cluster.vcov(x, survey3.2$state)))

model12.16<-lm(cmte_castetrust_n~treated+cmte_hh_index+contact_avg+
                 female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                 caste_reservation+state, survey3.2)
(coef12.16<-coeftest(model12.16, vcov=function(x) cluster.vcov(x, survey3.2$state)))

model12.18<-lm(neighbor_comfort_n~treated+cmte_hh_index+contact_avg+
                 female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                 caste_reservation+state, survey3.2)
(coef12.18<-coeftest(model12.18, vcov=function(x) cluster.vcov(x, survey3.2$state)))


model12.20<-lm(talk_comfort_n~treated+cmte_hh_index+contact_avg+
                 female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                 caste_reservation+state, survey3.2)
(coef12.20<-coeftest(model12.20, vcov=function(x) cluster.vcov(x, survey3.2$state)))

model12.22<-lm(donation_caste~treated+cmte_hh_index+contact_avg+
                 female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                 caste_reservation+state, survey3.2)
(coef12.22<-coeftest(model12.22, vcov=function(x) cluster.vcov(x, survey3.2$state)))


#Linear Models SI.6.5
stargazer(model12, model12.2, model12.4, model12.6, model12.8, model12.10, model12.12,
          model12.14, model12.16, model12.18, model12.20, model12.22,
          se=list(coef12[,2], coef12.2[,2], coef12.4[,2], coef12.6[,2], coef12.8[,2], coef12.10[,2],
                  coef12.12[,2], coef12.14[,2], coef12.16[,2], coef12.18[,2], coef12.20[,2], coef12.22[,2]))


#Table SI.6.3
#Forward Backward
model13<-lm(affect_mr_b2_p~treated+contactf_avg+cmte_pct_notforward+
              female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
              caste_reservation+corporation, survey3.2)
(coef13<-coeftest(model13, vcov=function(x) cluster.vcov(x, survey3.2$corporation)))


model13.2<-lm(affect_mr_b2_u~treated+contactf_avg+cmte_pct_notforward+
                female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                caste_reservation+corporation, survey3.2)
(coef13.2<-coeftest(model13.2, vcov=function(x) cluster.vcov(x, survey3.2$corporation)))

model13.4<-lm(affect_mr_b2_m~treated+contactf_avg+cmte_pct_notforward+
                female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                caste_reservation+corporation, survey3.2)
(coef13.4<-coeftest(model13.4, vcov=function(x) cluster.vcov(x, survey3.2$corporation)))

model13.6<-lm(affect_mr_b2_w~treated+contactf_avg+cmte_pct_notforward+
                female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                caste_reservation+corporation, survey3.2)
(coef13.6<-coeftest(model13.6, vcov=function(x) cluster.vcov(x, survey3.2$corporation)))

model13.8<-lm(cmte_spend_wisely_n~treated+contactf_avg+cmte_pct_notforward+
                female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                caste_reservation+corporation, survey3.2)
(coef13.8<-coeftest(model13.8, vcov=function(x) cluster.vcov(x, survey3.2$corporation)))

model13.10<-lm(cmte_opinions_n~treated+contactf_avg+cmte_pct_notforward+
                 female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                 caste_reservation+corporation, survey3.2)
(coef13.10<-coeftest(model13.10, vcov=function(x) cluster.vcov(x, survey3.2$corporation)))

model13.12<-lm(cmte_validconcerns_n~treated+contactf_avg+cmte_pct_notforward+
                 female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                 caste_reservation+corporation, survey3.2)
(coef13.12<-coeftest(model13.12, vcov=function(x) cluster.vcov(x, survey3.2$corporation)))

model13.14<-lm(cmte_onegroup_n~treated+contactf_avg+cmte_pct_notforward+
                 female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                 caste_reservation+corporation, survey3.2)
(coef13.14<-coeftest(model13, vcov=function(x) cluster.vcov(x, survey3.2$corporation)))

model13.16<-lm(cmte_castetrust_n~treated+contactf_avg+cmte_pct_notforward+
                 female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                 caste_reservation+corporation, survey3.2)
(coef13.16<-coeftest(model13.16, vcov=function(x) cluster.vcov(x, survey3.2$corporation)))


model13.18<-lm(neighbor_comfort_n~treated+contactf_avg+cmte_pct_notforward+
                 female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                 caste_reservation+corporation, survey3.2)
(coef13.18<-coeftest(model13.18, vcov=function(x) cluster.vcov(x, survey3.2$corporation)))

model13.20<-lm(talk_comfort_n~treated+contactf_avg+cmte_pct_notforward+
                 female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                 caste_reservation+corporation, survey3.2)
(coef13.20<-coeftest(model13.20, vcov=function(x) cluster.vcov(x, survey3.2$corporation)))

model13.22<-lm(donation_caste~treated+contactf_avg+cmte_pct_notforward+
                 female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                 caste_reservation+corporation, survey3.2)
(coef13.22<-coeftest(model13.22, vcov=function(x) cluster.vcov(x, survey3.2$corporation)))

#Linear Models SI.6.3
stargazer(model13, model13.2, model13.4, model13.6, model13.8, model13.10, model13.12,
          model13.14, model13.16, model13.18, model13.20, model13.22,
          se=list(coef13[,2], coef13.2[,2], coef13.4[,2], coef13.6[,2], coef13.8[,2], coef13.10[,2],
                  coef13.12[,2], coef13.14[,2], coef13.16[,2], coef13.18[,2], coef13.20[,2], coef13.22[,2]))



#Table SI.6.4
###Self-reported
model14<-lm(affect_mr_b2_p~treated+contact_avg_sr+cmte_hh_index_sr+
              female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
              caste_reservation+corporation, survey3.2)
(coef14<-coeftest(model14, vcov=function(x) cluster.vcov(x, survey3.2$corporation)))


model14.2<-lm(affect_mr_b2_u~treated+contact_avg_sr+cmte_hh_index_sr+
                female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                caste_reservation+corporation, survey3.2)
(coef14.2<-coeftest(model14.2, vcov=function(x) cluster.vcov(x, survey3.2$corporation)))

model14.4<-lm(affect_mr_b2_m~treated+contact_avg_sr+cmte_hh_index_sr+
                female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                caste_reservation+corporation, survey3.2)
(coef14.4<-coeftest(model14.4, vcov=function(x) cluster.vcov(x, survey3.2$corporation)))

model14.6<-lm(affect_mr_b2_w~treated+contact_avg_sr+cmte_hh_index_sr+
                female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                caste_reservation+corporation, survey3.2)
(coef14.6<-coeftest(model14.6, vcov=function(x) cluster.vcov(x, survey3.2$corporation)))

model14.8<-lm(cmte_spend_wisely_n~treated+contact_avg_sr+cmte_hh_index_sr+
                female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                caste_reservation+corporation, survey3.2)
(coef14.8<-coeftest(model14.8, vcov=function(x) cluster.vcov(x, survey3.2$corporation)))

model14.10<-lm(cmte_opinions_n~treated+contact_avg_sr+cmte_hh_index_sr+
                 female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                 caste_reservation+corporation, survey3.2)
(coef14.10<-coeftest(model14.10, vcov=function(x) cluster.vcov(x, survey3.2$corporation)))

model14.12<-lm(cmte_validconcerns_n~treated+contact_avg_sr+cmte_hh_index_sr+
                 female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                 caste_reservation+corporation, survey3.2)
(coef14.12<-coeftest(model14.12, vcov=function(x) cluster.vcov(x, survey3.2$corporation)))

model14.14<-lm(cmte_onegroup_n~treated+contact_avg_sr+cmte_hh_index_sr+
                 female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                 caste_reservation+corporation, survey3.2)
(coef14.14<-coeftest(model14, vcov=function(x) cluster.vcov(x, survey3.2$corporation)))

model14.16<-lm(cmte_castetrust_n~treated+contact_avg_sr+cmte_hh_index_sr+
                 female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                 caste_reservation+corporation, survey3.2)
(coef14.16<-coeftest(model14.16, vcov=function(x) cluster.vcov(x, survey3.2$corporation)))

model14.18<-lm(neighbor_comfort_n~treated+contact_avg_sr+cmte_hh_index_sr+
                 female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                 caste_reservation+corporation, survey3.2)
(coef14.18<-coeftest(model14.18, vcov=function(x) cluster.vcov(x, survey3.2$corporation)))

model14.20<-lm(talk_comfort_n~treated+contact_avg_sr+cmte_hh_index_sr+
                 female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                 caste_reservation+corporation, survey3.2)
(coef14.20<-coeftest(model14.20, vcov=function(x) cluster.vcov(x, survey3.2$corporation)))

model14.22<-lm(donation_caste~treated+contact_avg_sr+cmte_hh_index_sr+
                 female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                 caste_reservation+corporation, survey3.2)
(coef14.22<-coeftest(model14.22, vcov=function(x) cluster.vcov(x, survey3.2$corporation)))

#Linear Models SI.6.4
stargazer(model14, model14.2, model14.4, model14.6, model14.8, model14.10, model14.12,
          model14.14, model14.16, model14.18, model14.20, model14.22,
          se=list(coef14[,2], coef14.2[,2], coef14.4[,2], coef14.6[,2], coef14.8[,2], coef14.10[,2],
                  coef14.12[,2], coef14.14[,2], coef14.16[,2], coef14.18[,2], coef14.20[,2], coef14.22[,2]))



#Table SI.6.7
###Forward caste only
survey3.2$forward<-ifelse(survey3.2$Brahmin==1,1,NA)
survey3.2$forward<-ifelse(survey3.2$OF==1,1,survey3.2$forward)

survey3.2f<-survey3.2[survey3.2$forward==1,]

model15<-lm(affect_mr_b2_p~treated+contactf_avg+cmte_pct_notforward+
              female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
              corporation, survey3.2f)
(coef15<-coeftest(model15, vcov=function(x) cluster.vcov(x, survey3.2f$corporation)))

model15.2<-lm(affect_mr_b2_u~treated+contactf_avg+cmte_pct_notforward+
                female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                corporation, survey3.2f)
(coef15.2<-coeftest(model15.2, vcov=function(x) cluster.vcov(x, survey3.2f$corporation)))

model15.4<-lm(affect_mr_b2_m~treated+contactf_avg+cmte_pct_notforward+
                female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                corporation, survey3.2f)
(coef15.4<-coeftest(model15.4, vcov=function(x) cluster.vcov(x, survey3.2f$corporation)))

model15.6<-lm(affect_mr_b2_w~treated+contactf_avg+cmte_pct_notforward+
                female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                corporation, survey3.2f)
(coef15.6<-coeftest(model15.6, vcov=function(x) cluster.vcov(x, survey3.2f$corporation)))

model15.8<-lm(cmte_spend_wisely_n~treated+contactf_avg+cmte_pct_notforward+
                female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                corporation, survey3.2f)
(coef15.8<-coeftest(model15.8, vcov=function(x) cluster.vcov(x, survey3.2f$corporation)))

model15.10<-lm(cmte_opinions_n~treated+contactf_avg+cmte_pct_notforward+
                 female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                 corporation, survey3.2f)
(coef15.10<-coeftest(model15.10, vcov=function(x) cluster.vcov(x, survey3.2f$corporation)))

model15.12<-lm(cmte_validconcerns_n~treated+contactf_avg+cmte_pct_notforward+
                 female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                 corporation, survey3.2f)
(coef15.12<-coeftest(model15.12, vcov=function(x) cluster.vcov(x, survey3.2f$corporation)))

model15.14<-lm(cmte_onegroup_n~treated+contactf_avg+cmte_pct_notforward+
                 female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                 corporation, survey3.2f)
(coef15.14<-coeftest(model15, vcov=function(x) cluster.vcov(x, survey3.2f$corporation)))

model15.16<-lm(cmte_castetrust_n~treated+contactf_avg+cmte_pct_notforward+
                 female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                 corporation, survey3.2f)
(coef15.16<-coeftest(model15.16, vcov=function(x) cluster.vcov(x, survey3.2f$corporation)))

model15.18<-lm(neighbor_comfort_n~treated+contactf_avg+cmte_pct_notforward+
                 female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                 corporation, survey3.2f)
(coef15.18<-coeftest(model15.18, vcov=function(x) cluster.vcov(x, survey3.2f$corporation)))

model15.20<-lm(talk_comfort_n~treated+contactf_avg+cmte_pct_notforward+
                 female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                 corporation, survey3.2f)
(coef15.20<-coeftest(model15.20, vcov=function(x) cluster.vcov(x, survey3.2f$corporation)))

model15.22<-lm(donation_caste~treated+contactf_avg+cmte_pct_notforward+
                 female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                 corporation, survey3.2f)
(coef15.22<-coeftest(model15.22, vcov=function(x) cluster.vcov(x, survey3.2f$corporation)))


#Linear Models SI.6.7
stargazer(model15, model15.2, model15.4, model15.6, model15.8, model15.10, model15.12,
          model15.14, model15.16, model15.18, model15.20, model15.22,
          se=list(coef15[,2], coef15.2[,2], coef15.4[,2], coef15.6[,2], coef15.8[,2], coef15.10[,2],
                  coef15.12[,2], coef15.14[,2], coef15.16[,2], coef15.18[,2], coef15.20[,2], coef15.22[,2]))



#Table SI.6.6
###Now add SC/ST Violence
ncrb<-read.csv("NCRB_2017.csv", header=T, stringsAsFactors = F)
survey3.3<-merge(survey3.2, ncrb, by="corporation", all.x=T, all.y=F)
survey3.3$SC_ST_TotalCrime_log<-log(survey3.3$SC_ST_TotalCrime)

model18<-lm(affect_mr_b2_p~treated+cmte_hh_index+contact_avg+
              female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
              caste_reservation+SC_ST_TotalCrime_log+corporation, survey3.3)
(coef18<-coeftest(model18, vcov=function(x) cluster.vcov(x, survey3.3$corporation)))

model18.2<-lm(affect_mr_b2_u~treated+cmte_hh_index+contact_avg+
                female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                caste_reservation+SC_ST_TotalCrime_log+corporation, survey3.3)
(coef18.2<-coeftest(model18.2, vcov=function(x) cluster.vcov(x, survey3.3$corporation)))

model18.4<-lm(affect_mr_b2_m~treated+cmte_hh_index+contact_avg+
                female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                caste_reservation+SC_ST_TotalCrime_log+corporation, survey3.3)
(coef18.4<-coeftest(model18.4, vcov=function(x) cluster.vcov(x, survey3.3$corporation)))

model18.6<-lm(affect_mr_b2_w~treated+cmte_hh_index+contact_avg+
                female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                caste_reservation+SC_ST_TotalCrime_log+corporation, survey3.3)
(coef18.6<-coeftest(model18.6, vcov=function(x) cluster.vcov(x, survey3.3$corporation)))

model18.8<-lm(cmte_spend_wisely_n~treated+cmte_hh_index+contact_avg+
                female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                caste_reservation+SC_ST_TotalCrime_log+corporation, survey3.3)
(coef18.8<-coeftest(model18.8, vcov=function(x) cluster.vcov(x, survey3.3$corporation)))

model18.10<-lm(cmte_opinions_n~treated+cmte_hh_index+contact_avg+
                 female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                 caste_reservation+SC_ST_TotalCrime_log+corporation, survey3.3)
(coef18.10<-coeftest(model18.10, vcov=function(x) cluster.vcov(x, survey3.3$corporation)))

model18.12<-lm(cmte_validconcerns_n~treated+cmte_hh_index+contact_avg+
                 female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                 caste_reservation+SC_ST_TotalCrime_log+corporation, survey3.3)
(coef18.12<-coeftest(model18.12, vcov=function(x) cluster.vcov(x, survey3.3$corporation)))

model18.14<-lm(cmte_onegroup_n~treated+cmte_hh_index+contact_avg+
                 female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                 caste_reservation+SC_ST_TotalCrime_log+corporation, survey3.3)
(coef18.14<-coeftest(model18, vcov=function(x) cluster.vcov(x, survey3.3$corporation)))

model18.16<-lm(cmte_castetrust_n~treated+cmte_hh_index+contact_avg+
                 female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                 caste_reservation+SC_ST_TotalCrime_log+corporation, survey3.3)
(coef18.16<-coeftest(model18.16, vcov=function(x) cluster.vcov(x, survey3.3$corporation)))
marginal_temp<-as.data.frame(summary(margins::margins(model18.16, 
                                                      variables=c("treated"), 
                                                      vcov=vcovCL(model18.16, cluster=survey3.3$corporation))))
marginal_temp$dv<-"CasteTrust"
marginal<-rbind(marginal, marginal_temp)

model18.18<-lm(neighbor_comfort_n~treated+cmte_hh_index+contact_avg+
                 female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                 caste_reservation+SC_ST_TotalCrime_log+corporation, survey3.3)
(coef18.18<-coeftest(model18.18, vcov=function(x) cluster.vcov(x, survey3.3$corporation)))

model18.20<-lm(talk_comfort_n~treated+cmte_hh_index+contact_avg+
                 female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                 caste_reservation+SC_ST_TotalCrime_log+corporation, survey3.3)
(coef18.20<-coeftest(model18.20, vcov=function(x) cluster.vcov(x, survey3.3$corporation)))

model18.22<-lm(donation_caste~treated+cmte_hh_index+contact_avg+
                 female+age+socialmedia_active+edu_high+years_6over+BJP+INC+times_called_3over+
                 caste_reservation+SC_ST_TotalCrime_log+corporation, survey3.3)
(coef18.22<-coeftest(model18.22, vcov=function(x) cluster.vcov(x, survey3.3$corporation)))


#Linear Models SI.6.6
stargazer(model18, model18.2, model18.4, model18.6, model18.8, model18.10, model18.12,
          model18.14, model18.16, model18.18, model18.20, model18.22,
          se=list(coef18[,2], coef18.2[,2], coef18.4[,2], coef18.6[,2], coef18.8[,2], coef18.10[,2],
                  coef18.12[,2], coef18.14[,2], coef18.16[,2], coef18.18[,2], coef18.20[,2], coef18.22[,2]))






################################
###Main Text Plot: Figure 1
#load("marginal7.RData")

marginal_all<-marginal7
marginal_all$factor <- factor(marginal_all$factor, levels = c("contact_avg", "treated", "cmte_hh_index"),
                              labels = c("Contact", "Attention to Diversity", "Numeric Diversity"))
marginal_all$dv<-factor(marginal_all$dv,levels = c("Donation", "Talk", "Neighbor", "CasteTrust", "OneGroup",
                                                   "ValidConcerns", "Opinions", "SpendWisely", "Weak",
                                                   "Mixed", "Unpleasant", "Positive"),
                        labels=c("Donation", "Talk",
                                 "Neighbor", "Trust", "1 Group",
                                 "Valid", "Opinions", "Spend", "Weak",
                                 "Mixed", "Unpleas.", "Pleasant"))
marginal_all$type<-c(rep("Affect", 12), rep("Perceptions", 15), rep("Outgroup Attitudes", 9))
marginal_all$type<-factor(marginal_all$type, levels=c("Affect", "Perceptions", "Outgroup Attitudes"))


marginal_all1 <- ggplot(data=marginal_all, aes(dv, AME, color = factor, linetype=factor,
                                               group = factor, shape = factor)) + 
  geom_point(stat="identity", position = position_dodge(width = .4, 
                                                        preserve = "total"), size=3) + theme_bw()+
  facet_wrap(~type, strip.position = "top", scales = "free_y")


marginal_all2 <- marginal_all1 + geom_linerange(data=marginal_all, 
                                                aes(ymin=lower,ymax=upper), position = position_dodge(width = .4, preserve = "total"), size=1.3) + 
  labs(x="",y="Estimate") +
  scale_y_continuous(breaks=c(-0.6, -0.3, 0, 0.3, 0.6))+
  coord_flip(ylim=c(-0.65,0.65)) +  theme(axis.title=element_text(size=12)) + theme(axis.text=element_text(size=10)) + 
  theme(legend.text=element_text(size=12)) + theme(legend.title=element_text(size=12)) + 
  geom_hline(yintercept = 0, lty=2) +  
  theme(legend.position = "bottom") +theme(strip.text.x = element_text(size = 12))+
  labs(color  = "", linetype = "", shape = "")+
  scale_linetype_manual(values= c("solid", "solid", "solid"))+
  scale_colour_manual(values = c("orange3", "black", "gray45"))+
  scale_shape_manual(values = c(17,16,16))+
  theme(plot.margin = unit(c(0,0.3,0,-0.4), "cm"))
marginal_all2







